library(sfsmisc)

alpha<-.5;min<-1;max<-10000;n<-100000
setwd("C:\\Users\\Ofer\\Desktop\\My Dropbox\\PhD\\LATEX\\4\\figures")

f.dist<-function(x, min=1, max=10, alpha=1.4){alpha*x^(-alpha-1)/(min^-alpha - max^-alpha)}
f.cum<-function(x, min=1, max=10, alpha=1.4){(x^-alpha-max^-alpha)/(min^-alpha - max^-alpha)}
f.cum.app<-function(x, min=1, max=10, alpha=1.4){x^-alpha*(min^alpha*max^alpha)/(max^alpha - min^alpha)}
generatePower<- function(nmax, alpha, min=1, max=Inf) {
	u<-runif(nmax) # alpha<-1; min=1; max=10; u<-.001
	norm<-(min^-alpha - max^-alpha)
	(min^-alpha - u*norm)^-(1/alpha)
}
d<-generatePower(n, alpha, min=min, max=max)
# hist<-hist(d, plot = TRUE, breaks=n/10, freq=F, xlim=c(1,20), add=T, border="black")
hist.n<-hist(d, plot = FALSE, breaks=n,  freq=F, add=F, border="blue")
normalized<-rev(cumsum(rev(hist$density))) 


########################################################
x<-seq(1,12,length.out=100)
x.dist<-(hist.n$mid)[1:100]
y<-f.dist(x, min=min, max=max, alpha=alpha)
y.dist<-hist.n$density[1:100]
points((x), (y), pch=19, cex=.5)
points((x.dist), (y.dist), pch=19, col="blue",  cex=.2)

plot(log(x), log(y), pch=19, cex=.5)
points(log(x.dist), log(y.dist), pch=19, col="blue",  cex=.2)


lm<-lm(log(hist.n$density)[1:100] ~ log(hist.n$mid)[1:100])
lm<-lm(y.dist~log(x.dist))
summary(lm)
x<-seq(1,12,length.out=100)
y<-f.dist(x, min=min, max=max, alpha=alpha)
plot(log(x), log(y))
lm<-lm( log(y) ~ log(x))
summary(lm)



###########################################
# exact function of the distribution
plot(seq(min, max, length.out=n), f.dist(seq(min, max, length.out=n), min=min, max=max, alpha=alpha), 
	# from=min, to=max, 
	col="black", lwd=2,  type='l',
    	xlim=c(0,20),
	ylab='',  lty = 2, 
	xlab='', pch=20, cex=.1,  xaxt = "n",  yaxt = "n")

mtext("frequency", side=2, line=2.5, cex=1.5)
mtext("x-value", side=1, line=2.5, cex=1.5)
eaxis(1, cex.axis=.8)
eaxis(2, cex.axis=.8)

# integrate(f.dist, lower=min, upper=max, min=min, max=max, alpha=alpha)
# plot distribution
points( hist$mids, hist$density, 	col="blue", pch=3, cex=1, lwd=1.5)
#legend
leg.txt <- c("Exact function (DF)",
		 "Simulation     (DF)")
legend(x=15, y=.5, leg.txt,  lty = c(2, -1), pch = c(-1, 3),
      col = c("black", "blue"), 
	cex = 1, lwd=c(2,1))

dev.off()
